function [GlomStats] = AnalyseThresholdStatistics(stored_diams, stored_hist, PixelDimension, ThresholdRange, diams_range,MRorCT)
% Find the estimates for glomeruli number, size and distribution from
% threshold statistics of the depth image.  Return a struct
% containing :
%    optimal_index  index in the ThresholdRange with maximum num of gloms in first-pass
%    optimal_threshold threshold value of optimum glom count
%    GlomVoxelOptimum stuct for voxel method
%            mean    mean voxel size of clusters
%            sd      standard deviation of voxel size 
%            varinx  index to selected variable in fitobj 
%            total   sum of clusters within certain size range
%            histrange range of cluster sizes for counting total
%            fitobj  fitness result object 
%            gof     goodness-of-fit statistics
%            output  struct containing information associated with the fitting algorithm

%     Copyright © 2012-2013 Michael Eager <michael.eager@monash.edu>
%
%     This file is part of Xglom.
% 
%     This is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

if nargin < 6
    if PixelDimension> 0.04
        MRorCT=0;
    else
        MRorCT=1;
    end
        
end

if MRorCT == 0
    %% MR method
    %areas_estimate = sum(stored_hist(2:20,:),2);
    GlomStats.areas_estimate = sum(stored_hist(:,2:25)');
    Areas_minmax= extrems(GlomStats.areas_estimate);
    
    if  ~isempty(Areas_minmax.maxx)
        max_i = Areas_minmax.maxx(1);
        max_d = Areas_minmax.maxy(1);
    else
        [max_d max_i] = max(GlomStats.areas_estimate);
    end
    GlomStats.optimal_index = max_i
    GlomStats.optimal_threshold = ThresholdRange(max_i);
    
    
    opt_hist = stored_hist(max_i,1:100);
    
    x=1:numel(opt_hist);
    [fvoxels, gof, output] = fit(x(:),opt_hist(:),'gauss3')
    ci = confint(fvoxels);
    
    if (fvoxels.b2 < 20 && fvoxels.b2 > 3) && (ci(1,5) > 3 && ci(2,5) < 20 )
        GlomStats.GlomVoxelOptimum.mean = fvoxels.b2;
        GlomStats.GlomVoxelOptimum.sd = sqrt((fvoxels.c2)/2);
        GlomStats.GlomVoxelOptimum.varindx=5;
    elseif  (fvoxels.b1 < 20 && fvoxels.b1 > 3) && (ci(1,2) > 3 && ci(2,2) < 20 )
        GlomStats.GlomVoxelOptimum.mean = fvoxels.b1;
        GlomStats.GlomVoxelOptimum.sd = sqrt((fvoxels.c1)/2);
        GlomStats.GlomVoxelOptimum.varindx=2;
    elseif  (fvoxels.b3 < 20 && fvoxels.b3 > 3) && (ci(1,8) > 3 && ci(2,8) < 20 )
        GlomStats.GlomVoxelOptimum.mean = fvoxels.b3;
        GlomStats.GlomVoxelOptimum.sd = sqrt((fvoxels.c3)/2);
        GlomStats.GlomVoxelOptimum.varindx=8;
    else
        disp('Fitting procedure failed to find glom size range between 300 and 40 voxels')
        GlomStats.GlomVoxelOptimum.mean = 0;
        GlomStats.GlomVoxelOptimum.sd = 0;
        GlomStats.GlomVoxelOptimum.varindx=0;
    end
    
    % Estimate number of gloms from 2 voxels to 20 voxels
    GlomStats.GlomVoxelOptimum.total = sum(opt_hist(2:20));
    GlomStats.GlomVoxelOptimum.histrange = [2 20];
    
    GlomStats.GlomVoxelOptimum.fitobj = fvoxels;
    GlomStats.GlomVoxelOptimum.gof = gof;
    GlomStats.GlomVoxelOptimum.output=output;

    
    % Diamter estimation

    non_zero =stored_diams(max_i,:) >0;
    x=diams_range*PixelDimension*1000;%(non_zero)
    opt_diams = stored_diams(max_i,:);%non_zero)
    %opt_diams = smooth(opt_diams,5,'lowess');
    
    options = fitoptions('gauss3')
    options.Weights = non_zero;
    options.Lower = [0 0 0 0 0 0 0 0 0]
    %options.Robust = 'on'
    [fdiams,gof,output] = fit(x(:),opt_diams(:),'gauss3',options)
    ci=confint(fdiams);
    
    if  (fdiams.b1 < 200 && fdiams.b1 > 40) && (ci(1,2) > 40 && ci(2,2) < 200 )
        GlomStats.GlomDiamOptimum.mean = fdiams.b1;
        GlomStats.GlomDiamOptimum.sd = sqrt((fdiams.c1)/2);
        GlomStats.GlomDiamOptimum.varindx=2;
    elseif (fdiams.b2 < 200 && fdiams.b2 > 40) && (ci(1,5) > 40 && ci(2,5) < 200 )
        GlomStats.GlomDiamOptimum.mean = fdiams.b2;
        GlomStats.GlomDiamOptimum.sd = sqrt((fdiams.c2)/2);
        GlomStats.GlomDiamOptimum.varindx=5;
    elseif  (fdiams.b3 < 200 && fdiams.b3 > 40) && (ci(1,8) > 40 && ci(2,8) < 200 )
        GlomStats.GlomDiamOptimum.mean = fdiams.b3;
        GlomStats.GlomDiamOptimum.sd = sqrt((fdiams.c3)/2);
        GlomStats.GlomDiamOptimum.varindx=8;
    else
        disp('Fitting procedure failed to find glom diam size range between 40 and 200 microns')
        GlomStats.GlomDiamOptimum.mean = 0;
        GlomStats.GlomDiamOptimum.sd = 0;
        GlomStats.GlomDiamOptimum.varindx=0;
    end
    GlomStats.GlomDiamOptimum.total = sum(opt_diams(1:50));%(200:1000))
    GlomStats.GlomDiamOptimum.histrange = diams_range([1 50])*PixelDimension*1000;
    GlomStats.GlomDiamOptimum.fitobj = fdiams;
    GlomStats.GlomDiamOptimum.gof = gof;
    GlomStats.GlomDiamOptimum.output=output;
    

else
    %% micro CT
    % First pass -- Find the  best threshold assuming the gloms are around 100 voxels in size
    
    GlomStats.areas_estimate = sum(stored_hist(:,80:150)');
    
    Areas_minmax= extrems(GlomStats.areas_estimate);
    
    
    if  ~isempty(Areas_minmax.maxx)
        max_i = Areas_minmax.maxx(1);
        max_d = Areas_minmax.maxy(1);
    else
        [max_d max_i] = max(GlomStats.areas_estimate);
    end
    GlomStats.optimal_index = max_i
    GlomStats.optimal_threshold = ThresholdRange(max_i);
    
    % Find the mean and distribution
    % this process should be done using a fitting method
    %   opt_hist = stored_hist(max_i,:)
    %   hSmooth = smooth(opt_hist,10,'lowess',10);hSmooth = smooth(hSmooth,10);
    %   GlomAreas_minmax = extrems(hSmooth)
    %
    % if  ~isempty(GlomAreas_minmax.maxx) && ~isempty(GlomAreas_minmax.minx)
    %
    %
    %    GlomOptimum_mean = GlomAreas_minmax.maxx(1);
    %    GlomOptimim_sd = GlomAreas_minmax.maxy(1);
    %
    %    GlomOptimum_sum = stored_hist(max_i,GlomAreas_minmax.minx(1): 300)
    %
    % else
    %     [max_d max_i] = max(areas_estimate);
    % end
    
    opt_hist = stored_hist(max_i,1:end-1);
    opt_hist = smooth(opt_hist,10,'rlowess');
    %opt_hist = smooth(opt_hist,10,'rlowess');
    x=1:numel(opt_hist);
    options = fitoptions('gauss3');
    options.Lower = [0 -Inf 0 0 0 0 0 150 0];
    options.Upper = [Inf 0 Inf Inf 150 Inf Inf Inf Inf];
    [fvoxels,gof,output] = fit(x(:),opt_hist(:),'gauss3',options)
    ci=confint(fvoxels);
    
    g1 =  (fvoxels.a1).*exp(-((x-fvoxels.b1)./(fvoxels.c1)).^2);
    g2 =  (fvoxels.a2).*exp(-((x-fvoxels.b2)./(fvoxels.c2)).^2);
    g3 =  (fvoxels.a3).*exp(-((x-fvoxels.b3)./(fvoxels.c3)).^2);
    %hStat = figure; plot(x,g1,'-g',x,g2,'-b',x,g3,'-r', x,g1+g2+g3,'-k')
    
    if (fvoxels.b2 < 250 && fvoxels.b2 > 40) && (ci(1,5) > 40 && ci(2,5) < 250 )
        GlomStats.GlomVoxelOptimum.mean = fvoxels.b2;
        GlomStats.GlomVoxelOptimum.sd = sqrt((fvoxels.c2)/2);
        GlomStats.GlomVoxelOptimum.varindx=5;
    elseif  (fvoxels.b1 < 250 && fvoxels.b1 > 40) && (ci(1,2) > 40 && ci(2,2) < 250 )
        GlomStats.GlomVoxelOptimum.mean = fvoxels.b1;
        GlomStats.GlomVoxelOptimum.sd = sqrt((fvoxels.c1)/2);
        GlomStats.GlomVoxelOptimum.varindx=2;
    elseif  (fvoxels.b3 < 250 && fvoxels.b3 > 40) && (ci(1,8) > 40 && ci(2,8) < 250 )
        GlomStats.GlomVoxelOptimum.mean = fvoxels.b3;
        GlomStats.GlomVoxelOptimum.sd = sqrt((fvoxels.c3)/2);
        GlomStats.GlomVoxelOptimum.varindx=8;
    else
        disp('Fitting procedure failed to find glom size range between 300 and 40 voxels')
        GlomStats.GlomVoxelOptimum.mean = 0;
        GlomStats.GlomVoxelOptimum.sd = 0;
        GlomStats.GlomVoxelOptimum.varindx=0;
    end
    GlomStats.GlomVoxelOptimum.total = sum(opt_hist(24:250));
    GlomStats.GlomVoxelOptimum.histrange = [24 250];
    
    GlomStats.GlomVoxelOptimum.fitobj = fvoxels;
    GlomStats.GlomVoxelOptimum.gof = gof;
    GlomStats.GlomVoxelOptimum.output=output;
    
    
    
    
    
    x=diams_range(11:300)*PixelDimension*1000;
    non_zero =stored_diams(max_i,11:300) >0;
    opt_diams = stored_diams(max_i,11:300);
    opt_diams = smooth(opt_diams,10);%,'lowess');
    options = fitoptions('gauss3');
    options.Weights = non_zero;
    options.Lower = [0 0 0 0 0 0 0 0 0];
    [fdiams,gof,output] = fit(x(:),opt_diams(:),'gauss3',options)
    ci = confint(fdiams);
     GlomStats.GlomDiamOptimum.fitobj = fdiams;
    GlomStats.GlomVoxelOptimum.gof = gof;
    GlomStats.GlomVoxelOptimum.output=output;
    
    if  (fdiams.b3 < 200 && fdiams.b3 > 40) && (ci(1,8) > 40 && ci(2,8) < 200 )
        GlomStats.GlomDiamOptimum.mean = fdiams.b3;
        GlomStats.GlomDiamOptimum.sd = sqrt((fdiams.c3)/2);
        GlomStats.GlomDiamOptimum.varindx=8;
    elseif  (fdiams.b1 < 200 && fdiams.b1 > 40) && (ci(1,2) > 40 && ci(2,2) < 200 ) 
        GlomStats.GlomDiamOptimum.mean = fdiams.b1;
        GlomStats.GlomDiamOptimum.sd = sqrt((fdiams.c1)/2);
        GlomStats.GlomDiamOptimum.varindx=2;
    elseif (fdiams.b2 < 200 && fdiams.b2 > 40) && (ci(1,5) > 40 && ci(2,5) < 200 )
        GlomStats.GlomDiamOptimum.mean = fdiams.b2;
        GlomStats.GlomDiamOptimum.sd = sqrt((fdiams.c2)/2);
        GlomStats.GlomDiamOptimum.varindx=5;
    else
        disp('Fitting procedure failed to find glom diam size range between 40 and 200 microns')
        GlomStats.GlomDiamOptimum.mean = 0;
        GlomStats.GlomDiamOptimum.sd = 0;
        GlomStats.GlomDiamOptimum.varindx=0;
    end
    GlomStats.GlomDiamOptimum.total = sum(opt_diams(15:150));
    GlomStats.GlomDiamOptimum.histrange = [x(15) x(150)];
%    GlomStats.GlomDiamOptimum.histrange = diams_range([20 100])*PixelDimension*1000;


    
    
end
disp(GlomStats.GlomVoxelOptimum)
disp(GlomStats.GlomDiamOptimum)

% if handles.debug==1
%  plot(stored_hist(max_i,:)); xlim ([1 400])
%   hold('on'); %plot(hSmooth,'-r')
%   plot(GlomAreas_minmax.maxx,GlomAreas_minmax.maxy,'*r','Markersize',10)
%   plot(GlomAreas_minmax.minx,GlomAreas_minmax.miny,'*g','Markersize',10)
%   hold('off')
% end

